{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# DoWhy example on the Lalonde dataset\n",
    "\n",
    "Thanks to [@mizuy](https://github.com/mizuy) for providing this example. Here we use the Lalonde dataset and apply IPW estimator to it. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "R[write to console]: Loading required package: MASS\n",
      "\n",
      "R[write to console]: ## \n",
      "##  Matching (Version 4.9-7, Build Date: 2020-02-05)\n",
      "##  See http://sekhon.berkeley.edu/matching for additional documentation.\n",
      "##  Please cite software as:\n",
      "##   Jasjeet S. Sekhon. 2011. ``Multivariate and Propensity Score Matching\n",
      "##   Software with Automated Balance Optimization: The Matching package for R.''\n",
      "##   Journal of Statistical Software, 42(7): 1-52. \n",
      "##\n",
      "\n",
      "\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array(['Matching', 'MASS', 'tools', 'stats', 'graphics', 'grDevices',\n",
       "       'utils', 'datasets', 'methods', 'base'], dtype='<U9')"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import os, sys\n",
    "sys.path.append(os.path.abspath(\"../../../\"))\n",
    "\n",
    "import dowhy\n",
    "from dowhy import CausalModel\n",
    "from rpy2.robjects import r as R\n",
    "%load_ext rpy2.ipython\n",
    "\n",
    "#%R install.packages(\"Matching\")\n",
    "%R library(Matching)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Load the data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "%R data(lalonde)\n",
    "%R -o lalonde\n",
    "lalonde = lalonde.astype({'treat':'bool'}, copy=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Run DoWhy analysis: model, identify, estimate"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/amit/py-envs/env3.8/lib/python3.8/site-packages/sklearn/utils/validation.py:72: DataConversionWarning: A column-vector y was passed when a 1d array was expected. Please change the shape of y to (n_samples, ), for example using ravel().\n",
      "  return f(**kwargs)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Causal Estimate is 1639.8075601428254\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<caption>WLS Regression Results</caption>\n",
       "<tr>\n",
       "  <th>Dep. Variable:</th>          <td>re78</td>       <th>  R-squared:         </th> <td>   0.015</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Model:</th>                   <td>WLS</td>       <th>  Adj. R-squared:    </th> <td>   0.013</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Method:</th>             <td>Least Squares</td>  <th>  F-statistic:       </th> <td>   6.743</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Date:</th>             <td>Sun, 06 Jun 2021</td> <th>  Prob (F-statistic):</th>  <td>0.00973</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Time:</th>                 <td>19:11:54</td>     <th>  Log-Likelihood:    </th> <td> -4544.7</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>No. Observations:</th>      <td>   445</td>      <th>  AIC:               </th> <td>   9093.</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Df Residuals:</th>          <td>   443</td>      <th>  BIC:               </th> <td>   9102.</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Df Model:</th>              <td>     1</td>      <th>                     </th>     <td> </td>   \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Covariance Type:</th>      <td>nonrobust</td>    <th>                     </th>     <td> </td>   \n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "        <td></td>           <th>coef</th>     <th>std err</th>      <th>t</th>      <th>P>|t|</th>  <th>[0.025</th>    <th>0.975]</th>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Intercept</th>     <td> 4555.0759</td> <td>  406.704</td> <td>   11.200</td> <td> 0.000</td> <td> 3755.767</td> <td> 5354.385</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>treat[T.True]</th> <td> 1639.8076</td> <td>  631.496</td> <td>    2.597</td> <td> 0.010</td> <td>  398.708</td> <td> 2880.907</td>\n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "  <th>Omnibus:</th>       <td>303.262</td> <th>  Durbin-Watson:     </th> <td>   2.085</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Prob(Omnibus):</th> <td> 0.000</td>  <th>  Jarque-Bera (JB):  </th> <td>4770.581</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Skew:</th>          <td> 2.709</td>  <th>  Prob(JB):          </th> <td>    0.00</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Kurtosis:</th>      <td>18.097</td>  <th>  Cond. No.          </th> <td>    2.47</td>\n",
       "</tr>\n",
       "</table><br/><br/>Notes:<br/>[1] Standard Errors assume that the covariance matrix of the errors is correctly specified."
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.summary.Summary'>\n",
       "\"\"\"\n",
       "                            WLS Regression Results                            \n",
       "==============================================================================\n",
       "Dep. Variable:                   re78   R-squared:                       0.015\n",
       "Model:                            WLS   Adj. R-squared:                  0.013\n",
       "Method:                 Least Squares   F-statistic:                     6.743\n",
       "Date:                Sun, 06 Jun 2021   Prob (F-statistic):            0.00973\n",
       "Time:                        19:11:54   Log-Likelihood:                -4544.7\n",
       "No. Observations:                 445   AIC:                             9093.\n",
       "Df Residuals:                     443   BIC:                             9102.\n",
       "Df Model:                           1                                         \n",
       "Covariance Type:            nonrobust                                         \n",
       "=================================================================================\n",
       "                    coef    std err          t      P>|t|      [0.025      0.975]\n",
       "---------------------------------------------------------------------------------\n",
       "Intercept      4555.0759    406.704     11.200      0.000    3755.767    5354.385\n",
       "treat[T.True]  1639.8076    631.496      2.597      0.010     398.708    2880.907\n",
       "==============================================================================\n",
       "Omnibus:                      303.262   Durbin-Watson:                   2.085\n",
       "Prob(Omnibus):                  0.000   Jarque-Bera (JB):             4770.581\n",
       "Skew:                           2.709   Prob(JB):                         0.00\n",
       "Kurtosis:                      18.097   Cond. No.                         2.47\n",
       "==============================================================================\n",
       "\n",
       "Notes:\n",
       "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n",
       "\"\"\""
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model=CausalModel(\n",
    "        data = lalonde,\n",
    "        treatment='treat',\n",
    "        outcome='re78',\n",
    "        common_causes='nodegr+black+hisp+age+educ+married'.split('+'))\n",
    "identified_estimand = model.identify_effect(proceed_when_unidentifiable=True)\n",
    "estimate = model.estimate_effect(identified_estimand,\n",
    "        method_name=\"backdoor.propensity_score_weighting\",\n",
    "        target_units=\"ate\",                         \n",
    "        method_params={\"weighting_scheme\":\"ips_weight\"})\n",
    "#print(estimate)\n",
    "print(\"Causal Estimate is \" + str(estimate.value))\n",
    "\n",
    "import statsmodels.formula.api as smf\n",
    "reg=smf.wls('re78~1+treat', data=lalonde, weights=lalonde.ips_stabilized_weight)\n",
    "res=reg.fit()\n",
    "res.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Interpret the estimate\n",
    "The plot below shows how the distribution of a confounder, \"married\" changes from the original data to the weighted data. In both datasets, we compare the distribution of \"married\" across treated and untreated units."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAsgAAAHwCAYAAAC7apkrAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAyKklEQVR4nO3de5hcZZXv8e8iAYIkco2cmBCDGjXgwSAhIEYmXAYQVFABwTMQECeCUQcRFRlHooNHdFSO4BFkhCGoKAhyiIoXrgIiYNCIQGCIEiAhQAwQQLmFrPPH3p286VR3upOurr58P8/TT1ftveuttXdVr/717reqIjORJEmSVNmg1QVIkiRJfYkBWZIkSSoYkCVJkqSCAVmSJEkqGJAlSZKkggFZkiRJKgzYgBwR50TEv/XQWGMj4pmIGFJfvz4iPtgTY9fj/TwipvXUeN2439Mi4q8R8Uhv33cjEfG2iLh3HW87NSIWdrI+I+K16zj2WyPivvo5cPC6jNEXdPV5uz6Pw7pYn8emGSLiroiY2sVtF0TEPs2taOCxP3fpfgdNf17LbSMi/isinoiI29ZljFay73ZNX+y7/TIg1wfn2Yh4OiKejIibI+K4iFi5P5l5XGb+exfH6vRAZ+aDmTk8M1/qgdpnRsT32o3/9syctb5jd7OOscAngO0z83/05n13JDNvzMzXt7qOBr4AfLN+Dvy/VhfTbM18HHo6vDRDZu6Qmdev7zjrEwr6M/vz+rM/r2YK8I/AmMyc3OgxGgjsu32v7/bLgFx7Z2aOAF4FnA58Gjivp+8kIob29Jh9xFhgaWY+1tt33OiY9vHj/CrgrnW5YR/fL6lZ7M/rx/68yquABZn5t54YbAA/Z9TTMrPffQELgH3aLZsMrADeWF+/ADitvrw18FPgSeBx4EaqPw6+W9/mWeAZ4FPAOCCBY4EHgRuKZUPr8a4HvgTcBjwFXAFsWa+bCixsVC+wP/AC8GJ9f38sxvtgfXkD4LPAA8BjwIXAZvW6tjqm1bX9FfjXTo7TZvXtl9TjfbYef596n1fUdVzQ4LZTgYX1MXkMWAwcDBwA/Hd9HE9pd/x/Wx/jxcA3gY2K9QnMAO4D7i/G/zTwSP1YrHbsgFcCl9X13w98rFi3Sf0YPwHcDXyy/XFvtz8JfAz4S33c/gPYoFj/AWBePd4vgVfVy//M6s+Rjeu6ZtfHYD7wz8U4M4FLge/Vz40P1o/DefVxWQScBgzpoM6uHMeG+wEcDfymvs0y4B5g7+K219f1bFTX/j+Lda8A/g6MbPA4LABOAu6ox70YGFas/1Rd68P1+Am8tsG+fRF4CXiuPpbfLPbpuPq58STwf4FY22PTYPxZwCfqy6PrcWfU119T73PbsXoHMLe+v5uBHRv1F6rn2az6vufV+7rWYwNsyuo/Y89QPW8mA3Pq58ajwNdb3U97+gv7s/25+/35G8BD9eN1O/C2evmxVP3ipfpY3NrBY9Rhj2VVXzwDWEr9vGvw/LTv2ndX37dWN9N1+aJBA66XPwgcX1++gFUN+EvAOcCG9dfb2p4I7cdiVZO7sD7Ym9C4AS8C3lhvcxnwvbJxdVQvVYD6Xrv117OqAX+AKnS9GhgO/Bj4brva/rOu603A88CEDo7ThVS/HEbUt/1v4NiO6mx326nAcuBz9TH7Z6pGeFE93g71E3G7evudgd2AofV9zQNOaNdgrgK2rGtvG//LVKGzbdnCevsNqBrl56gay6upmtN+9frTqX6RbglsC9y5lv1J4Lp6+7H1sWg75gfVx3xCXf9ngZs7er5R/VL+FtUP5MT6uOxVPL4vUv2y2qDer8uBb1M9V15B9Yv7Qx3U2ZXj2NF+HF0f04/Xj9n7qJpHWzi4vtj2W8CXi3H/BfhJo+dGvf+3UTWaLeuajqvX7U/1C3QH4GVUfxg0bNTta2i3Tz8FNq/3aQmwf1cem3bjfKDYh/dT/XFzcbHuivryTlShYldgCFWgWQBs3ODn9XTg18AWwBiqhtzVY7PacayX/RY4sr48HNit1f20p7+wP9ufu9+f/wnYqq7vE1Q9ZVi97mjgpmLbRo9Rhz2WVX3xo/X4mzS4f/uufXfNfevt5tkTX3TcgG+h/oud1RvwF6gaUaO/rlYbi1VN7tUNlpUN+PRi/fZUf9UO6eDBKR/4mXTegK8BPlysez1V4Gr7wU2quVht628DDm+wX0PqmrYvln0IuL6jJ1G720+larBtf4WPqO9712Kb24GDO7j9CcDl7X4Y92o3/gus/hfxypqofogebDfmZ4D/qi//hfqHub4+fS37k+22/zBwTX3559S/mOrrG1D9Vf+qBo/ftlR/jY8otv8S9Vme+vG9oVi3DdUvyU2KZUcA13Xxud7oOHa0H0dTnU0ozwLcxqrGUD7PdqUKLG1BZA5wWKPnRr3//1Rc/wpwTn35fOBLxbrXsm6Nekpx/RLg5K48Nu3GeQ3VGYcNqALXh4rn0yzgxPry2cC/t7vtvcA/NHi8V/7Sr69/sBvHZrXjWC+7Afg8sHVXHv/++IX92f7czf7coL4ngDfVl4+mk4DMWnpsffsHu3rfnRwf++4g67v9eQ5yI6OpTue39x9Ufw39KiL+EhEnd2Gsh7qx/gGqvxy37lKVnXtlPV459lCqJtCmfFXz36n+Impv67qm9mON7kYtS3PVC1+erb8/Wqx/tu2+I+J1EfHTiHgkIp4C/jdrHo/2x3RJZj7XwX2/Cnhl/SKfJyPiSeAUVh2HV7LmY7A27bd/ZXFf3yju53EgaHysXgk8nplPtxur3La8n1dRPQ6Li/G/TXWWYw3rcBzL/QBYlHVH6GA9AJl5K9VzZ2pEvIGqwc5uVFOto+dc+8dhbT833R2/y49NZv4Z+BvVWf23UZ0deTgiXg/8A9UZibYxP9HuubUtDY4TXdu/rvw8tjkWeB1wT0T8LiLe0cm2A439eRX7cyEiToqIeRGxrB5rswb1daQrPbbT54t9d43x7bv07xfprSYidqF68G5qvy4zn87MT2Tmq4F3ASdGxN5tqzsYsqPlbbYtLo+lOovwV6onysuKuoZQzS/q6rgPUz2RyrGXs3rj64q/1jW1H2tRN8fpqrOp5l6Nz8yXUzXLaLdN+33v7Fg8BNyfmZsXXyMy84B6/WLWfAzWpv32Dxf39aF297VJZt7cYIyHgS0jYkS7scrjWu7XQ1RnN7Yuxn55Zu7QQY1dOY4d7QfA6IiITtaXZlH9a/NI4NJOfhl2ZjHVv8Aa1dbI2p7/7XXnsYGqGR9CNX9wUX19GtW/6uYWY36x3Zgvy8wfNBivu/tXWmNfM/O+zDyC6pf3l4FLI2LTbozZL9mf12B/rkXE26jmmB4GbJGZm1NNUWhfX0d1daXHru1xte+uzr7LAAjIEfHy+q+BH1L92+VPDbZ5R0S8tn4CL6P6F/mKevWjVPOnuuufImL7iHgZ1b8IL63/mv9vYFhEHBgRG1LN3dm4uN2jwLjyLY/a+QHw8YjYLiKGU/0le3FmLu9OcXUtlwBfjIgREfEq4ESquUrNMIJqAvwz9V/Gx6/neLcBT0fEpyNik4gYEhFvrH/RQrVvn4mILSJiDNX8srX5ZL39tlRzvy6ul59Tj7UDQERsFhGHNhogMx+ienHBlyJiWETsSPXXacPjmpmLgV8BX6ufqxtExGsi4h86qLErx7Gj/YCqAXwsIjas92ECcGUH9/U94N1UzfrCDrZZm0uAYyJiQv2zsLb3tu3uz1uXH5var4GPUP1LDap/LX6E6l+0bWfb/hM4LiJ2jcqm9c/riDWHW+15Nroeq6seBbaKiM3aFkTEP0XEyMxcQfVCFVjViwYc+3Nj9uc1altONQd2aER8Dnh5J9uv9hitQ4/tqAb77ir2Xfp3QP5JRDxN9VfJvwJfB47pYNvxwNVUr2j8LfCtzLyuXvcl4LNRnfI/qRv3/12qeXSPUL1Y62MAmbmMan7Sd6jOBvyN6tXAbX5Uf18aEb9vMO759dg3UL0y+Dm6Fv4a+Wh9/3+hOnNzUT1+M5xENUH/aaofhIs737xz9Q/VO6j+bXM/1RmX71D96w2q+UQP1Ot+RXXM1uYKqnl5c4GfUb/tVGZeTvVX5Q+j+vfancDbOxnnCKr5hg9TvTjk1My8upPtj6J6IcvdVHO1LgVGdbBtV45jw/2o3Ur1fP8r1auXD8nMpY3uqA77v6f6i/vGTurvUGb+HDiT6gUs86nmmUJ1RqeRbwCHRPWm/2d2YfzuPja/pvpl19aob6I6Y9h2ncycQ/Wipm9SPR7zqeYRNvIFqp/f+6l6yKWd7Fv72u+hClR/qfvLK6leXHNXRDxDdSwOz8xnOxunn7I/r539ufJL4BdUf7w8QHVMO5sy0Ogx6k6PbcS+u/r49l1WTRSX1A9ERFL9G3B+g3VHU70QY0o3xjsfeDgzP9tD9U2gaqYbd/esWn8QEcdTNdfunJ2S1I/Zd1urVX23P59BlrQeImIc8B7W8wMcIuLdEbFxRGxBddbhJwOlSUfEqKg+anyDqF508gmq/xpIUrfZd9eur/RdA7I0CEXEv1OdcfiPzLx/PYf7ENX7W/6Zav7o+s5v7Es2onpF/NPAtVT/Zv1WSyuS1C/Zd7usT/Rdp1hIkiRJBc8gS5IkSYWhrS5gfWy99dY5bty4VpchSd12++23/zUzR659y77PXiypv+qoF/frgDxu3DjmzJnT6jIkqdsioiuf/tgv2Isl9Vcd9WKnWEiSJEkFA7IkSZJUMCBLkiRJhX49B1lSz3rxxRdZuHAhzz33XKtLGTCGDRvGmDFj2HDDDVtdiqR+wD7cHN3txQZkSSstXLiQESNGMG7cOCKi1eX0e5nJ0qVLWbhwIdttt12ry5HUD9iHe9669GKnWEha6bnnnmOrrbayKfeQiGCrrbbyTJCkLrMP97x16cUGZEmrsSn3LI+npO6yb/S87h5TA7IkSZJUcA6ypA6NO/lnPTregtMPXPs2Cxbwjne8gzvvvHPlspkzZzJ8+HBOOumkhreZO3cuDz/8MAcccECP1Lmu402dOpWvfvWrTJo0qUfqkCT7cGv6sGeQJfV7c+fO5corr2y4bvny5T06niRpTQOtDxuQJfUbU6dO5dOf/jSTJ0/mda97HTfeeCMvvPACn/vc57j44ouZOHEiF198MTNnzuTII4/krW99K0ceeSRLlizhve99L7vssgu77LILv/nNbwC47bbbeMtb3sJOO+3E7rvvzr333ttwvL/97W984AMfYPLkyey0005cccUVADz77LMcfvjhTJgwgXe/+908++yzrTw8ktR0g6UPO8VCUr+yfPlybrvtNq688ko+//nPc/XVV/OFL3yBOXPm8M1vfhOo/hV49913c9NNN7HJJpvw/ve/n49//ONMmTKFBx98kP3224958+bxhje8gRtvvJGhQ4dy9dVXc8opp3DZZZetMd4pp5zCXnvtxfnnn8+TTz7J5MmT2Wefffj2t7/Ny172MubNm8cdd9zBm9/85lYeGknqFYOhDxuQJfUpHb3SuG35e97zHgB23nlnFixY0OE473rXu9hkk00AuPrqq7n77rtXrnvqqad45plnWLZsGdOmTeO+++4jInjxxRcbjvWrX/2K2bNn89WvfhWo3obpwQcf5IYbbuBjH/sYADvuuCM77rhj93ZWkvog+7ABWVIfs9VWW/HEE0+stuzxxx9f+ebuG2+8MQBDhgzpdF7bpptuuvLyihUruOWWWxg2bNhq23zkIx9hzz335PLLL2fBggVMnTq14ViZyWWXXcbrX//6ddklSepX7MPOQZbUxwwfPpxRo0Zx7bXXAlVT/sUvfsGUKVM6vM2IESN4+umnO1y/7777ctZZZ628PnfuXACWLVvG6NGjAbjgggs6HG+//fbjrLPOIjMB+MMf/gDAHnvswUUXXQTAnXfeyR133NGNPZWkvsk+7BlkSZ3oytsBNcOFF17IjBkzOPHEEwE49dRTec1rXtPh9nvuuSenn346EydO5DOf+cwa688880xmzJjBjjvuyPLly9ljjz0455xz+NSnPsW0adM47bTTOPDAAzsc79/+7d844YQT2HHHHVmxYgXbbbcdP/3pTzn++OM55phjmDBhAhMmTGDnnXfu+YMhaVCzD7emD0dbEu+PJk2alHPmzGl1GdKAMW/ePCZMmNDqMgacRsc1Im7PzAHxhsn2Yqnn2Iebpzu92CkWkiRJUsGALEmSJBUG5Rzknv7Yxr6oVXOWJKmrBnovtg9L/ZdnkCVJkqSCAVmSJEkqDMopFpIkNd3MzVpdQfPNXNbqCqSmMCBL6lhP/4Jfyy/TpUuXsvfeewPwyCOPMGTIEEaOHAnAbbfdxkYbbdTtu7z++uvZaKON2H333bt1u3HjxjFnzhy23nrrbt+nJPUY+3BL+rABWVKfsdVWW638dKWZM2cyfPhwTjrppJXrly9fztCh3Wtb119/PcOHD+92Y5akwcg+XHEOsqQ+7eijj+a4445j11135VOf+hR//vOf2X///dl5551529vexj333APAT37yE3bddVd22mkn9tlnHx599FEWLFjAOeecwxlnnMHEiRO58cYbWbJkCe9973vZZZdd2GWXXfjNb34DVGdN9t13X3bYYQc++MEP0p8/REmSetJg7MOeQZbU5y1cuJCbb76ZIUOGsPfee3POOecwfvx4br31Vj784Q9z7bXXMmXKFG655RYigu985zt85Stf4Wtf+xrHHXfcamdA3v/+9/Pxj3+cKVOm8OCDD7Lffvsxb948Pv/5zzNlyhQ+97nP8bOf/YzzzjuvxXstSX3HYOvDBmRJfd6hhx7KkCFDeOaZZ7j55ps59NBDV657/vnngap5v+9972Px4sW88MILbLfddg3Huvrqq7n77rtXXn/qqad45plnuOGGG/jxj38MwIEHHsgWW2zRxD2SpP5lsPVhA7KkPm/TTTcFYMWKFWy++eYr58eVPvrRj3LiiSfyrne9i+uvv56ZM2c2HGvFihXccsstDBs2rIkVS9LAMtj6sHOQJfUbL3/5y9luu+340Y9+BEBm8sc//hGAZcuWMXr0aABmzZq18jYjRozg6aefXnl933335ayzzlp5va3J77HHHlx00UUA/PznP+eJJ55o6r5IUn80WPqwZ5AldawPvsfp97//fY4//nhOO+00XnzxRQ4//HDe9KY3MXPmTA499FC22GIL9tprL+6//34A3vnOd3LIIYdwxRVXcNZZZ3HmmWcyY8YMdtxxR5YvX84ee+zBOeecw6mnnsoRRxzBDjvswO67787YsWNbvKeShH24RaI/v1J70qRJOWfOnG7fbtzJP2tCNX3LgtMPbHUJ6ofmzZvHhAkTWl3GgNPouEbE7Zk5qUUl9Sh7cWMLhr2/1SU0Xx8Mb/2dfbh5utOLnWIhSZIkFQzIkiRJUsGALGk1/XnaVV/k8ZTUXfaNntfdY2pAlrTSsGHDWLp0qc25h2QmS5cu7dNvZSSpb7EP97x16cW+i4WklcaMGcPChQtZsmRJq0sZMIYNG8aYMWNaXYakfsI+3Bzd7cUGZEkrbbjhhh1+8pH6p4gYBtwAbEzV8y/NzFMj4gLgH4C2tyE4OjPnRkQA3wAOAP5eL/9971cuDU724b7BgCxJA9vzwF6Z+UxEbAjcFBE/r9d9MjMvbbf924Hx9deuwNn1d0kaNJyDLEkDWFaeqa9uWH91NrnxIODC+na3AJtHxKhm1ylJfYkBWZIGuIgYEhFzgceAqzLz1nrVFyPijog4IyI2rpeNBh4qbr6wXtZ+zOkRMSci5jhXUtJAY0CWpAEuM1/KzInAGGByRLwR+AzwBmAXYEvg090c89zMnJSZk0aOHNnTJUtSSxmQJWmQyMwngeuA/TNzcT2N4nngv4DJ9WaLgG2Lm42pl0nSoGFAlqQBLCJGRsTm9eVNgH8E7mmbV1y/a8XBwJ31TWYDR0VlN2BZZi7u9cIlqYV8FwtJGthGAbMiYgjVSZFLMvOnEXFtRIwEApgLHFdvfyXVW7zNp3qbt2N6v2RJai0DsiQNYJl5B7BTg+V7dbB9AjOaXZck9WVOsZAkSZIKBmRJkiSpYECWJEmSCgZkSZIkqWBAliRJkgoGZEmSJKlgQJYkSZIKBmRJkiSp0LSAHBHbRsR1EXF3RNwVEf9SL98yIq6KiPvq71vUyyMizoyI+RFxR0S8uVm1SZIkSR1p5hnk5cAnMnN7YDdgRkRsD5wMXJOZ44Fr6usAbwfG11/TgbObWJskSZLUUNMCcmYuzszf15efBuYBo4GDgFn1ZrOAg+vLBwEXZuUWYPOIGNWs+iRJkqRGemUOckSMA3YCbgW2yczF9apHgG3qy6OBh4qbLayXtR9rekTMiYg5S5YsaV7RkiRJGpSaHpAjYjhwGXBCZj5VrsvMBLI742XmuZk5KTMnjRw5sgcrlSRJkpockCNiQ6pw/P3M/HG9+NG2qRP198fq5YuAbYubj6mXSZIkSb2mme9iEcB5wLzM/HqxajYwrb48DbiiWH5U/W4WuwHLiqkYkiRJUq8Y2sSx3wocCfwpIubWy04BTgcuiYhjgQeAw+p1VwIHAPOBvwPHNLE2SZIkqaGmBeTMvAmIDlbv3WD7BGY0qx5JkiSpK/wkPUmSJKlgQJYkSZIKBmRJkiSpYECWJEmSCgZkSZIkqWBAliRJkgoGZEmSJKlgQJYkSZIKBmRJkiSpYECWJEmSCgZkSZIkqWBAliRJkgoGZEmSJKlgQJYkSZIKBmRJkiSpYECWJEmSCgZkSZIkqWBAliRJkgoGZEmSJKlgQJYkSZIKBmRJkiSpYECWJEmSCgZkSZIkqWBAliRJkgoGZEmSJKlgQJYkSZIKBmRJkiSpYECWJEmSCgZkSZIkqWBAliRJkgoGZEkawCJiWETcFhF/jIi7IuLz9fLtIuLWiJgfERdHxEb18o3r6/Pr9eNaugOS1AIGZEka2J4H9srMNwETgf0jYjfgy8AZmfla4Ang2Hr7Y4En6uVn1NtJ0qBiQJakASwrz9RXN6y/EtgLuLRePgs4uL58UH2dev3eERG9U60k9Q0GZEka4CJiSETMBR4DrgL+DDyZmcvrTRYCo+vLo4GHAOr1y4CtGow5PSLmRMScJUuWNHkPJKl3GZAlaYDLzJcycyIwBpgMvKEHxjw3Mydl5qSRI0eu73CS1KcYkCVpkMjMJ4HrgLcAm0fE0HrVGGBRfXkRsC1AvX4zYGnvVipJrWVAlqQBLCJGRsTm9eVNgH8E5lEF5UPqzaYBV9SXZ9fXqddfm5nZawVLUh8wdO2bSJL6sVHArIgYQnVS5JLM/GlE3A38MCJOA/4AnFdvfx7w3YiYDzwOHN6KoiWplQzIkjSAZeYdwE4Nlv+Faj5y++XPAYf2QmmS1Gc5xUKSJEkqGJAlSZKkggFZkiRJKhiQJUmSpIIBWZIkSSoYkCVJkqSCAVmSJEkqGJAlSZKkggFZkiRJKhiQJUmSpIIBWZIkSSoYkCVJkqSCAVmSJEkqGJAlSZKkggFZkiRJKhiQJUmSpIIBWZIkSSoYkCVJkqSCAVmSJEkqGJAlSZKkggFZkiRJKhiQJUmSpIIBWZIkSSoYkCVJkqSCAVmSJEkqGJAlSZKkggFZkiRJKhiQJUmSpIIBWZIkSSoYkCVJkqSCAVmSJEkqGJAlSZKkggFZkiRJKhiQJUmSpIIBWZIkSSoYkCVJkqSCAVmSJEkqGJAlSZKkggFZkiRJKhiQJUmSpIIBWZIkSSoYkCVJkqSCAVmSJEkqGJAlSZKkggFZkiRJKhiQJUmSpELTAnJEnB8Rj0XEncWymRGxKCLm1l8HFOs+ExHzI+LeiNivWXVJkiRJnWnmGeQLgP0bLD8jMyfWX1cCRMT2wOHADvVtvhURQ5pYmyQNChGxbURcFxF3R8RdEfEv9XJPWEhSB4Y2a+DMvCEixnVx84OAH2bm88D9ETEfmAz8tln1SdIgsRz4RGb+PiJGALdHxFX1ujMy86vlxu1OWLwSuDoiXpeZL/Vq1ZLUQk0LyJ34SEQcBcyhatpPAKOBW4ptFtbL1hAR04HpAGPHjm1yqf3YzM1aXUHzzVzW6gqkPi8zFwOL68tPR8Q8OuivNU9YSBr0evtFemcDrwEmUjXsr3V3gMw8NzMnZeakkSNH9nB5kjRw1f/V2wm4tV70kYi4o37NyBb1stHAQ8XNOjxhIUkDVa8G5Mx8NDNfyswVwH9SnZUAWARsW2w6pl4mSeoBETEcuAw4ITOfYj1PWETE9IiYExFzlixZ0tPlSlJL9WpAjohRxdV3A23vcDEbODwiNo6I7YDxwG29WZskDVQRsSFVOP5+Zv4Y1v+Ehf/NkzSQNW0OckT8AJgKbB0RC4FTgakRMRFIYAHwIYDMvCsiLgHupnpByQxfECJJ6y8iAjgPmJeZXy+Wj6rnJ8OaJywuioivU71IzxMWkgadZr6LxRENFp/XyfZfBL7YrHokaZB6K3Ak8KeImFsvOwU4whMWktRYK97FQpLUSzLzJiAarLqyk9t4wkLSoOZHTUuSJEkFA7IkSZJUMCBLkiRJBQOyJEmSVDAgS5IkSQUDsiRJklQwIEuSJEkFA7IkSZJUMCBLkiRJBQOyJEmSVDAgS5IkSQUDsiRJklQwIEuSJEkFA7IkSZJUMCBLkiRJBQOyJEmSVDAgS5IkSQUDsiRJklQwIEuSJEkFA7IkSZJUMCBLkiRJBQOyJEmSVDAgS5IkSQUDsiRJklQwIEuSJEkFA7IkSZJUMCBLkiRJBQOyJEmSVDAgS5IkSQUDsiRJklQwIEuSJEkFA7IkSZJUMCBLkiRJBQOyJEmSVDAgS5IkSQUDsiRJklQwIEuSJEkFA7IkSZJUMCBLkiRJBQOyJEmSVDAgS5IkSQUDsiRJklQwIEuSJEkFA7IkSZJUMCBLkiRJBQOyJEmSVOhSQI6It3ZlmSSpOezDktR7unoG+awuLpMkNYd9WJJ6ydDOVkbEW4DdgZERcWKx6uXAkGYWJkla/z4cEdsCFwLbAAmcm5nfiIgtgYuBccAC4LDMfCIiAvgGcADwd+DozPx9z+2RJPV9azuDvBEwnCpIjyi+ngIOaW5pkiTWvw8vBz6RmdsDuwEzImJ74GTgmswcD1xTXwd4OzC+/poOnN1zuyJJ/UOnZ5Az89fAryPigsx8oJdqkiTV1rcPZ+ZiYHF9+emImAeMBg4CptabzQKuBz5dL78wMxO4JSI2j4hR9TiSNCh0GpALG0fEuVT/ilt5m8zcqxlFSZLWsN59OCLGATsBtwLbFKH3EaopGFCF54eKmy2sl60WkCNiOtUZZsaOHduN3ZCkvq+rAflHwDnAd4CXmleOJKkD69WHI2I4cBlwQmY+VU01rmRmRkR2Z7zMPBc4F2DSpEnduq0k9XVdDcjLM9N5aJLUOuvchyNiQ6pw/P3M/HG9+NG2qRMRMQp4rF6+CNi2uPmYepkkDRpdfZu3n0TEhyNiVERs2fbV1MokSaV16sP1u1KcB8zLzK8Xq2YD0+rL04AriuVHRWU3YJnzjyUNNl09g9zWRD9ZLEvg1T1bjiSpA+vah98KHAn8KSLm1stOAU4HLomIY4EHgMPqdVdSvcXbfKq3eTtmvSuXpH6mSwE5M7drdiGSpI6tax/OzJuA6GD13g22T2DGutyXJA0UXQrIEXFUo+WZeWHPliNJasQ+LEm9p6tTLHYpLg+jOuvwe6pPZ5IkNZ99WJJ6SVenWHy0vB4RmwM/bEZBkqQ12Yclqfd09V0s2vsb4LxkSWod+7AkNUlX5yD/hOrV0gBDgAnAJc0qSpK0OvuwJPWers5B/mpxeTnwQGYubEI9kqTG7MOS1Eu6NMUiM38N3AOMALYAXmhmUZKk1dmHJan3dCkgR8RhwG3AoVRvJn9rRBzSzMIkSavYhyWp93R1isW/Artk5mMAETESuBq4tFmFSZJWYx+WpF7S1Xex2KCtKdeWduO2kqT1Zx+WpF7S1TPIv4iIXwI/qK+/D7iyOSVJkhqwD0tSL+k0IEfEa4FtMvOTEfEeYEq96rfA95tdnCQNdvZhSep9azuD/H+AzwBk5o+BHwNExP+s172zibVJkuzDktTr1jZ/bZvM/FP7hfWycU2pSJJUsg9LUi9bW0DevJN1m/RgHZKkxjbvZJ19WJKaYG0BeU5E/HP7hRHxQeD25pQkSSrYhyWpl61tDvIJwOUR8b9Y1YgnARsB725iXZKkygnYhyWpV3UakDPzUWD3iNgTeGO9+GeZeW3TK5Mk2YclqQW69D7ImXkdcF2Ta5EkdcA+LEm9x09hkiRJkgoGZEmSJKlgQJYkSZIKBmRJkiSpYECWJEmSCgZkSZIkqWBAliRJkgoGZEmSJKnQtIAcEedHxGMRcWexbMuIuCoi7qu/b1Evj4g4MyLmR8QdEfHmZtUlSZIkdaaZZ5AvAPZvt+xk4JrMHA9cU18HeDswvv6aDpzdxLokSZKkDjUtIGfmDcDj7RYfBMyqL88CDi6WX5iVW4DNI2JUs2qTJEmSOtLbc5C3yczF9eVHgG3qy6OBh4rtFtbL1hAR0yNiTkTMWbJkSfMqlSRJ0qDUshfpZWYCuQ63OzczJ2XmpJEjRzahMkmSJA1mvR2QH22bOlF/f6xevgjYtthuTL1MkiRJ6lW9HZBnA9Pqy9OAK4rlR9XvZrEbsKyYiiFJkiT1mqHNGjgifgBMBbaOiIXAqcDpwCURcSzwAHBYvfmVwAHAfODvwDHNqkuSJEnqTNMCcmYe0cGqvRtsm8CMZtUiSZIkdZWfpCdJkiQVDMiSJElSwYAsSZIkFQzIkiRJUsGALEmSJBUMyJIkSVLBgCxJkiQVDMiSJElSwYAsSZIkFQzIkiRJUsGALEmSJBUMyJI0gEXE+RHxWETcWSybGRGLImJu/XVAse4zETE/Iu6NiP1aU7UktZYBWZIGtguA/RssPyMzJ9ZfVwJExPbA4cAO9W2+FRFDeq1SSeojDMiSNIBl5g3A413c/CDgh5n5fGbeD8wHJjetOEnqowzIkjQ4fSQi7qinYGxRLxsNPFRss7BetoaImB4RcyJizpIlS5pdqyT1KgOyJA0+ZwOvASYCi4GvdXeAzDw3Mydl5qSRI0f2cHmS1FoGZEkaZDLz0cx8KTNXAP/JqmkUi4Bti03H1MskaVAxIEvSIBMRo4qr7wba3uFiNnB4RGwcEdsB44Hbers+SWq1oa0uQJLUPBHxA2AqsHVELAROBaZGxEQggQXAhwAy866IuAS4G1gOzMjMl1pQtiS1lAFZkgawzDyiweLzOtn+i8AXm1eRJPV9TrGQJEmSCgZkSZIkqWBAliRJkgoGZEmSJKlgQJYkSZIKBmRJkiSpYECWJEmSCgZkSZIkqWBAliRJkgoGZEmSJKlgQJYkSZIKBmRJkiSpYECWJEmSCgZkSZIkqWBAliRJkgoGZEmSJKlgQJYkSZIKBmRJkiSpYECWJEmSCgZkSZIkqWBAliRJkgoGZEmSJKlgQJYkSZIKBmRJkiSpYECWJEmSCgZkSZIkqWBAliRJkgoGZEmSJKlgQJYkSZIKBmRJkiSpYECWJEmSCgZkSZIkqWBAliRJkgoGZEmSJKlgQJYkSZIKBmRJkiSpYECWJEmSCgZkSZIkqWBAliRJkgoGZEmSJKlgQJYkSZIKBmRJkiSpYECWJEmSCgZkSZIkqWBAliRJkgoGZEmSJKlgQJYkSZIKBmRJGsAi4vyIeCwi7iyWbRkRV0XEffX3LerlERFnRsT8iLgjIt7cusolqXUMyJI0sF0A7N9u2cnANZk5Hrimvg7wdmB8/TUdOLuXapSkPsWALEkDWGbeADzebvFBwKz68izg4GL5hVm5Bdg8Ikb1SqGS1IcYkCVp8NkmMxfXlx8BtqkvjwYeKrZbWC9bQ0RMj4g5ETFnyZIlzatUklrAgCxJg1hmJpDrcLtzM3NSZk4aOXJkEyqTpNYxIEvS4PNo29SJ+vtj9fJFwLbFdmPqZZI0qBiQJWnwmQ1Mqy9PA64olh9Vv5vFbsCyYiqGJA0aQ1tdgCSpeSLiB8BUYOuIWAicCpwOXBIRxwIPAIfVm18JHADMB/4OHNPrBUtSH2BAlqQBLDOP6GDV3g22TWBGcyuSpL7PKRaSJElSwYAsSZIkFQzIkiRJUsGALEmSJBUMyJIkSVLBgCxJkiQVDMiSJElSwYAsSZIkFQzIkiRJUqEln6QXEQuAp4GXgOWZOSkitgQuBsYBC4DDMvOJVtQnSZKkwauVZ5D3zMyJmTmpvn4ycE1mjgeuqa9LkiRJvaovTbE4CJhVX54FHNy6UiRJkjRYtSogJ/CriLg9IqbXy7bJzMX15UeAbVpTmiRJkgazlsxBBqZk5qKIeAVwVUTcU67MzIyIbHTDOlBPBxg7dmzzK5UkSdKg0pIzyJm5qP7+GHA5MBl4NCJGAdTfH+vgtudm5qTMnDRy5MjeKlmSJEmDRK8H5IjYNCJGtF0G9gXuBGYD0+rNpgFX9HZtkiRJUiumWGwDXB4Rbfd/UWb+IiJ+B1wSEccCDwCHtaA2SZIkDXK9HpAz8y/AmxosXwrs3dv1SJIkSaW+9DZvkiRJUssZkCVJkqRCq97mTZIkSb1h5matrqD5Zi7r0eE8gyxJkiQVPIMsSZIGrXEn/6zVJTTdgmGtrqD/8QyyJEmSVDAgS5IkSQUDsiRJklQwIEuSJEkFA7IkSZJUMCBLkiRJBQOyJEmSVDAgS5IkSQU/KERqkUHx5vSnH9jqEiRJ6jbPIEuSJEkFA7IkSZJUMCBLkiRJBQOyJEmSVDAgS5IkSQUDsiRJklQwIEuSJEkFA7IkSZJUMCBLkiRJBT9JT1LzzNys1RU038xlra5AktTDPIMsSZIkFQzIkiRJUsGALEmSJBUMyJIkSVLBgCxJkiQVfBcLSRqkImIB8DTwErA8MydFxJbAxcA4YAFwWGY+0aoaJakVPIMsSYPbnpk5MTMn1ddPBq7JzPHANfV1SRpUDMiSpNJBwKz68izg4NaVIkmtYUCWpMErgV9FxO0RMb1etk1mLq4vPwJs05rSJKl1nIMsSYPXlMxcFBGvAK6KiHvKlZmZEZGNblgH6ukAY8eObX6lktSLPIMsSYNUZi6qvz8GXA5MBh6NiFEA9ffHOrjtuZk5KTMnjRw5srdKlqReYUCWpEEoIjaNiBFtl4F9gTuB2cC0erNpwBWtqVCSWscpFpI0OG0DXB4RUP0uuCgzfxERvwMuiYhjgQeAw1pYoyS1hAFZkgahzPwL8KYGy5cCe/d+RZLUdzjFQpIkSSoYkCVJkqSCAVmSJEkqGJAlSZKkggFZkiRJKhiQJUmSpIIBWZIkSSoYkCVJkqSCAVmSJEkqGJAlSZKkggFZkiRJKhiQJUmSpIIBWZIkSSoYkCVJkqSCAVmSJEkqGJAlSZKkggFZkiRJKhiQJUmSpIIBWZIkSSoYkCVJkqSCAVmSJEkqGJAlSZKkggFZkiRJKhiQJUmSpIIBWZIkSSoYkCVJkqSCAVmSJEkqGJAlSZKkggFZkiRJKhiQJUmSpIIBWZIkSSoYkCVJkqSCAVmSJEkqGJAlSZKkggFZkiRJKhiQJUmSpIIBWZIkSSoYkCVJkqSCAVmSJEkqGJAlSZKkggFZkiRJKhiQJUmSpIIBWZIkSSoYkCVJkqSCAVmSJEkqGJAlSZKkggFZkiRJKhiQJUmSpEKfC8gRsX9E3BsR8yPi5FbXI0mDjX1Y0mDXpwJyRAwB/i/wdmB74IiI2L61VUnS4GEflqQ+FpCBycD8zPxLZr4A/BA4qMU1SdJgYh+WNOgNbXUB7YwGHiquLwR2LTeIiOnA9PrqMxFxby/V1q8EbA38tdV1NNXno9UVaC18HnbqVT1ZRg9aax8Ge3FX+PxXX+FzsVMNe3FfC8hrlZnnAue2uo6+LiLmZOakVtehwc3n4cBlL147n//qK3wudl9fm2KxCNi2uD6mXiZJ6h32YUmDXl8LyL8DxkfEdhGxEXA4MLvFNUnSYGIfljTo9akpFpm5PCI+AvwSGAKcn5l3tbis/sp/faov8HnYz9iHe5TPf/UVPhe7KTKz1TVIkiRJfUZfm2IhSZIktZQBWZIkSSoYkPu5tX0kbERsHBEX1+tvjYhxLShTA1xEnB8Rj0XEnR2sj4g4s34e3hERb+7tGqVmsQ+rL7AP9ywDcj/WxY+EPRZ4IjNfC5wBfLl3q9QgcQGwfyfr3w6Mr7+mA2f3Qk1S09mH1YdcgH24xxiQ+7eufCTsQcCs+vKlwN4R4UcfqUdl5g3A451schBwYVZuATaPiFG9U53UVPZh9Qn24Z5lQO7fGn0k7OiOtsnM5cAyYKteqU5apSvPVak/sg+rv7APd4MBWZIkSSoYkPu3rnwk7MptImIosBmwtFeqk1bx44s1UNmH1V/Yh7vBgNy/deUjYWcD0+rLhwDXpp8Oo943GziqfhX1bsCyzFzc6qKkHmAfVn9hH+6GPvVR0+qejj4SNiK+AMzJzNnAecB3I2I+1eT9w1tXsQaqiPgBMBXYOiIWAqcCGwJk5jnAlcABwHzg78AxralU6ln2YfUV9uGe5UdNS5IkSQWnWEiSJEkFA7IkSZJUMCBLkiRJBQOyJEmSVDAgS5IkSQUDstRFEfHKiLi0m7e5ICIOaVZNkjTY2IvVGwzIUgP1p12tdj0zH85MG6wk9RJ7sVrFgKwBJSLGRcQ99dmC/46I70fEPhHxm4i4LyIm11+/jYg/RMTNEfH6+rZHR8TsiLgWuKbB9XERcWe97ZCI+I+I+F1E3BERH6qXR0R8MyLujYirgVe07GBIUovYi9Xf+Ul6GoheCxwKfIDqY2DfD0wB3gWcAhwFvK3+BKx9gP8NvLe+7ZuBHTPz8Yg4ut31ccV9HEv1MZ27RMTGwG8i4lfATsDrge2BbYC7gfObubOS1EfZi9VvGZA1EN2fmX8CiIi7gGsyMyPiT8A4YDNgVkSMB5L6ozhrV2Xm451cb7MvsGMxp20zYDywB/CDzHwJeLg+4yFJg5G9WP2WAVkD0fPF5RXF9RVUz/l/B67LzHfXZyKuL7b/W7ux2l9vE8BHM/OXqy2MOGAda5akgcZerH7LOcgajDYDFtWXj17HMX4JHB8RGwJExOsiYlPgBuB99by4UcCe61usJA1Q9mL1WQZkDUZfAb4UEX9g3f+L8h2qOW2/r18s8u16rMuB++p1FwK/Xf9yJWlAsherz4rMbHUNkiRJUp/hGWRJkiSpYECWJEmSCgZkSZIkqWBAliRJkgoGZEmSJKlgQJYkSZIKBmRJkiSp8P8BjfbB6AKQoTAAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 720x504 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "estimate.interpret(method_name=\"confounder_distribution_interpreter\",var_type='discrete',\n",
    "                   var_name='married', fig_size = (10, 7), font_size = 12)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Sanity check: compare to manual IPW estimate"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Causal Estimate is 1639.8075601428245\n"
     ]
    }
   ],
   "source": [
    "df = model._data\n",
    "ps = df['ps']\n",
    "y = df['re78']\n",
    "z = df['treat']\n",
    "\n",
    "ey1 = z*y/ps / sum(z/ps)\n",
    "ey0 = (1-z)*y/(1-ps) / sum((1-z)/(1-ps))\n",
    "ate = ey1.sum()-ey0.sum()\n",
    "print(\"Causal Estimate is \" + str(ate))\n",
    "\n",
    "# correct -> Causal Estimate is 1634.9868359746906"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.5"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": false,
   "sideBar": true,
   "skip_h1_title": true,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
